3D-SLIP model based dynamic stability strategy for legged robots with impact disturbance rejection

Inspired by biomechanical studies, the spring-loaded inverted pendulum model is an effective behavior model to describe the running movement of animals and legged robots in the sagittal plane. However, when confronted with external lateral disturbances, the model has to move out of the 2-D plane and be extended to 3-D locomotion. With the degree of freedom increasing, the computational complexity is higher and the real-time control is more and more difficult, especially when considering the complex legged model. Here, we construct a control strategy based on the classical Raibert controller for legged locomotion under lateral impact disturbances. This strategy, named 3D-HFC, is composed of three core modules: touchdown angle control, body attitude angle control and energy compensation. The energy loss in each step is taken into consideration, and the real-time measured energy loss of the current step is adopted to predict that of the next step. We demonstrate the efficiency of the proposed control strategy on a simulated 3D-SLIP lower order model and a simulated running quadruped, which are perturbed by different impact forces. Furthermore, a quadruped bionic prototype named MBBOT was set up, on which lateral impact experiments were designed and implemented. Both simulation and experimental results show that the proposed approach can realize the impact disturbance rejection.

It has been well-recognized that legged locomotion has superior performance over its wheeled and tracked counterparts on rough environment due to its discrete footholds, active suspension configuration and larger working space 1 . The past few decades have witnessed a large number of endeavors which were conducted to improve the dynamic performance of legged robots. The spring-loaded inverted pendulum model (SLIP) has been established as a template to approximate the dynamics of trotting, running, and hopping in 2-, 4-, 6-, and 8-legged animals 2,3 . One of the most representative examples employing this idea is the series of monopod and quadruped running machines built by Raibert and his colleagues 4 . After that, many other hopping robots have been built [5][6][7][8][9][10] , and some legged robots have been designed with the SLIP-like feature [11][12][13][14][15] .
In order to meet the requirements of the real-time control system, different online control methods have been proposed for the SLIP model in the sagittal plane. Raibert provided the first control algorithm for SLIP-based running robots 16 . Koditschek and Buhler introduced two models constrained to hop in one dimension using a simplified Raibert scheme 17 . Taking it one step further, Zeglin and Brown built a planar bow leg monopod robot and controlled it with a combination of graph-searching planning and feedback control algorithm 18 . Saranli et al. proposed control algorithms to make multi-jointed planar monopod hoppers behave like spring-mass systems 19 . Liu presented a control approach for bounding gait of quadruped robots with multi-jointed legs by applying the concept of Virtual Constraints 20 . Nie proposed an open-loop control method including touchdown angle control strategy and leg length control strategy that can drive the quadruped robot with spring legs running in passive bounding gait 21 . Buehler presented the CPDR controller imposing desired trajectories via inverse dynamics to reduce locomotion energy expenditure 22 . The CPDR strategy is able to realize the control of hopping height through energy regulation, but the main concerns of this adjustment are the robot's vertical velocity components. Yu developed a dead-beat controller for the SLIP runner to regulate the apex state based on approximate apex return map [23][24][25] . He also proposed Sagittal SLIP-anchored task space control for a monopode robot traversing irregular terrain 26 . Schmitt and Clark proposed the AER control strategy based on active energy addition and removal by leg compression and leg extension 27 . Arslan and Saranli realized dead-beat control of hopping machines using three control methods, i.e. the Leg Length Control (LLC), Leg Stiffness Control (LSC) and Two-Phase Stiffness Control (TPSC) 28 . This approach regulated system energy by changing leg length and Methods Equivalent 3D-SLIP model. Bionic principle. Multiped mammals have very complex leg structure, for example, a dog's hind leg is composed of femur, tibia, fibula, and hind paw, with muscles, tendons and ligaments providing traction forces for the dog's body. Complicated connection between bones, redundant muscle actuation, and excessive DOFs of joints all lead to high dimension, high computational complexity and high redundancy features of legged locomotion, making its analysis very difficult. Based on observation and analysis of various multipedes running at different speeds, biologists find in legged locomotion the property of elastic energy storage 41,42 , and therefore propose the Spring Loaded Inverted Pendulum (SLIP) model virtually composed of a mass and a leg spring to simplify the description of legged locomotion. The single leg simplification procedure is shown in Fig. 1.
Locomotion gaits of multiped mammals can be divided into symmetrical gaits and asymmetrical gaits 43 . Raibert put forward the monopod SLIP model in sagittal plane as the equivalence to rapid symmetric running gait 16 . Take dog's trotting gait as an example, its diagonal legs have consistent motion phase, touching down on and taking off from the ground at the same time, thus this pair of legs can be equivalently replaced by a single virtual leg at body CoM. The stiffness, force, and torque of the single virtual leg are the synthesis of those of the two diagonal legs. In trotting gait, two pairs of diagonal legs alternate in motion with phase difference to be one half, and the kinetic and dynamic characteristics of these two pairs are symmetric. Therefore, the monopod SLIP model equivalently placed at the body CoM can describe the motion characteristics of quadruped mammals. The establishment of the equivalent virtual leg is shown in Fig. 2 www.nature.com/scientificreports/ 3D-SLIP model. Control of the sagittal 2D-SLIP model generally involves only the control of one rotational DOF of the leg spring about the body, and the regulation of body sagittal pitch angle. To extend this model to three-dimensional space, another rotational DOF between body and leg spring is added, and three spatial attitude angles are taken into consideration 44 . The 3D-SLIP model proposed in this paper is shown as Fig. 3.
The 3D-SLIP model consists of a body with lumped mass and equivalent leg spring. in Fig. 3, let point O and point P denote the equivalent body CoM and the leg toe, respectively. Let e b be the fixed coordinate system defined at the body's CoM, and e b x , e b y , e b z as its base vector. Body and leg spring are connected with a two-DOF rotational hinge at equivalent body CoM. Link 1 has the rotational DOF around base vector e b x , with rotation angle defined as δ 1 . Link 2 has the rotational DOF in plane L 1 OL 2 with rotation angle defined as δ 2 . Let e r be the reference coordinate defined at leg toe, and e r x , e r y , e r z as its base vector. Body attitude is described with generalized Euler angles (α, β, γ ) , the sequence of which is a rotation of an angle α about the e b z axis, followed by a rotation β about the e b y axis, and then followed by a rotation γ about the e b x axis. BP, OA and AB are the projections of leg spring OP on base vectors e r z , e r y , e r x respectively. AOB is defined as side-swing plane, BOP is defined as forward-swing plane. Leg spring attitude is described with side and forward touchdown angles in reference coordinate, namely ϕ in AOB and θ in BOP. Equivalent body mass is defined as    www.nature.com/scientificreports/ m body , and the equivalent leg spring is considered massless to facilitate subsequent dynamic analysis and calculation. Leg length, equivalent stiffness, and equivalent damping are described with r, k and c. The SLIP model is actuated with three inputs: two torques τ 1 and τ 2 applied at the hip around z-axis and x-axis respectively and a force F acting along the leg.
Dynamic analysis of the 3D-SLIP model. For SLIP model running in sagittal plane, when confronted with a sudden lateral impact force, a velocity perpendicular to this motion plane will be produced, extending the locomotion of SLIP from 2D plane to 3D space. As shown in Fig. 4, in each step, the locomotion of SLIP model in three-dimensional space consists of flight phase and stance phase, with ϕ TD and θ TD denoting side and forward touchdown angles in TOUCHDOWN event respectively, ϕ LO and θ LO denoting side and forward liftoff angles in LIFTOFF event respectively. In flight phase, the APEX is triggered when the SLIP reaches its maximum height, where the vertical velocity decreases to zero. States of SLIP can be represented by the characteristic parameters of APEX event in each step, including the horizontal velocity ẋ , lateral velocity ż , vertical height of body CoM y APEX and system energy E. The apex return map is defined from the current apex variable U n [ẋ n ,ż n , y APEX n , E n ] to the next variable U n+1 [ẋ n+1 ,ż n+1 , y APEX n+1 , E n+1 ] , where n represents the number of steps counting from 1. There are two different dynamic characteristics for the 3D-SLIP model in flight phase and stance phase respectively. Because the system is limited only by the effect of gravity in flight phase, the flight dynamics of the model in Cartesian coordinates can be easily obtained by In stance phase, 3D-SLIP can be seen as an inverted pendulum with leg toe fixed at the ground. We define a polar coordinate system at equivalent body CoM along the direction of leg spring and the direction perpendicular to the leg plane. Body translational velocity can be decomposed into along the direction of leg spring, in leg plane and in side-swing plane perpendicular to the leg. The dynamics equation of stance phase can be derived as where r 0 is the original leg length, r is the real-time leg length in stance phase.
It is obvious that Eq. (2) is nonlinear differential equations with no exact analytical solutions. In 2D-SLIP cases, various methods have been proposed to search for approximate solutions like the linearization of gravity or using small enough compression ratio of leg spring 23,45,46 . Similar methods can be extended to 3D-SLIP cases, but to obtain solutions with reasonable precision for control with the employment of these methods requires more numerical computation, and will encounter more difficulties in the real-time control. Besides the environmental factors may also affect the real-time control of the system. But the environmental factors can be considered as the external disturbance uniformly. In this paper, a real-time control strategy is developed to realize a 3-D space periodical stable motion control under external impact forces.
Control method under the impact of disturbance. In previous study, Raibert first proposed a three-part control system to achieve the dynamically stable running for the hopping and quadruped running machines 4 . Based on the Raibert's controller, we have proposed the Hybrid Feedback Control (HFC) strategy for www.nature.com/scientificreports/ stability control of 2D-SLIP in sagittal plane on uneven terrain 30 . Further, by extending HFC strategy to threedimensional space and taking into consideration the control of body attitude angles, this paper puts forward a Three-Dimension Hybrid Feedback Control (3D-HFC) strategy for 3D-SLIP model 44 .
Control system structure for the 3D-HFC method. Figure 5 shows the control system structure of the 3D-HFC. It is composed of three core modules: touchdown angle control module (TAC), body attitude angle control module (AAC) and energy compensation control module (ECC). Choose U n [ẋ n ,ż n , y APEX n , E n ] as the APEX state variable characterizing the state of SLIP system in each step, U des [ẋ des ,ż des , y APEX des , E des ] denoting the desired state. The goal of the controller is to get U n as close as possible to the desired state U des . The current state U n can be measured through the sensor module at the every apex. Compared with the desired state, the relative error of the state can be obtained and the hybrid feedback control of the SLIP system can be completed through the above-mentioned three control modules.
Touchdown angle control strategy. In stance phase, 2D-SLIP can be seen as an inverted pendulum model. With initial conditions determined, this model can realize the periodical jumping from touchdown to liftoff through passive compression and release of the leg spring. 3D space locomotion is similar. Different from the 2D case, we have to control touchdown angles in two directions for 3D-SLIP, θ TD n in forward-swing plane and ϕ TD n in sideswing plane. 3D-SLIP touchdown angle control is shown in Fig. 6.
The projections of the touchdown point P in forward-swing plane and side-swing plane are defined as P TD x and P TD z respectively. The neutral points which have been introduced in Raibert's book 4 are defined as P N x and P N z in forward-swing plane and side-swing plane. S nx and S nz represent distances between projections of touchdown point and neutral point in two planes respectively. According to Raibert's foot placement algorithms 4 , S nx and S nz can be described as follows:  www.nature.com/scientificreports/ where µ 1 and µ 2 are feedback gains selected to maximize stability. S 0 nx and S 0 nz respectively correspond to the relative position of the neutral points P N x and P N z in two projection planes, and they are half the distances of the horizontal movement of the CoM during the stance phase on the condition it has achieved symmetry movement in two planes. S 0 nx and S 0 nz can be obtained through the product of the average speed and the stance phase time, and then they can be determined respectively as where ẋ avg n−1 and ż avg n−1 are the average forward speed and the average lateral speed during the stance phase, and T S is the approximate stance phase time 8 . r TD n−1 , θ TD n−1 , ϕ TD n−1 , t TD n−1 and r LO n−1 , θ LO n−1 , ϕ LO n−1 , t LO n−1 represent actual leg length, forward-swing angle, side-swing angle and time at TOUCHDOWN and LIFTOFF phase of the previous stance, respectively. They all can be obtained through the sensing system. Then the forward touchdown angle θ TD n and side touchdown angle ϕ TD n in Step n can be calculated as follows: Control strategy for body attitude angle. In the simplified SLIP model, the mass of the leg is ignored. With no external forces acting on the system, the body posture remains unchanged during its motion procedure. But in a real platform, the leg mass cannot be set to zero. This leads to a cumulative deviation of the body attitude. Besides, encountered with instantaneous lateral impact forces, body changes its posture. Therefore, the body attitude needs to be regulated in the overall stability control strategy.
Since system angular momentum is conservative in flight phase, the body attitude control could only be realized through the regulation of two driving torques at the hip in stance phase. For model simplification, we assume no slip between the leg toe and the ground in stance phase. The torques controlled to act on the hip around z-axis and x-axis respectively in this phase are determined with linear equations as follows: where k α,p and k γ ,p are position feedback gains, k α,v and k γ ,v are velocity feedback gains. We can determine the approximate range of these gains based on the simulation results, in order to ensure our first attempt will not damage the system. And then, we can optimize these gains through empirical tuning.
Energy compensation control strategy. Previous studies ignore the damping of leg spring in SLIP model and the energy loss in contact between leg toe and the ground. Some researchers take a step further to consider energy loss of SLIP running on flat terrain. However, adding a constant energy to the system for compensation is inapplicable for energy control of system running on rough terrain or under external impact forces. In order to reflect the moving process of practical running systems more precisely, the energy loss in each step is taken into consideration and the real-time measured energy loss of the current step is adopted to predict that of the next step, avoiding the complexity of precise computation. In the Step (n), E + n . and E − n denote the amount of energy added and lost during the stance phase respectively. After the liftoff event, the amount of energy loss E − n can be calculated as follows: where E TD n and E LO n are the overall system mechanical energy at the moment of touchdown and liftoff in the Step (n) respectively. They are determined as follows: where ẋ TD n ,ẏ TD n ,ż TD n , r TD n , θ TD n , ϕ TD n and ẋ LO n ,ẏ LO n ,ż LO n , r LO n , θ LO n , ϕ LO n represent the horizontal velocities, vertical velocities, lateral velocities, leg spring lengths, forward-swing touchdown/liftoff angles and side-swing touchdown/liftoff angles at the touchdown and liftoff moment respectively. These system state variables are obtained through the sensor module in simulation and the Inertial Measurement Units (IMU) on the robot platform.
When the system state converges to the steady-state, the fluctuation of the energy loss gets smaller. Therefore, it assumes that the energy loss of the steady system to be the same as the previous cycle. The amount of compensational energy in the Step (n + 1) can be calculated as www.nature.com/scientificreports/ Since the system is energy conservative in flight phase, the system energy should be compensated in the stance phase. Some scholars have proposed different energy compensation strategies like Rated Impulse control (RIC), Constant Force Control (CFC), Variable Stiffness Control (VSC) 23,24 . The VSC method is chosen in this paper because its mechanism is closer to that of practical system. The compensation force of VSC can be determined as where r B is the minimum spring length, F VSC n+1 acts on the body CoM from bottom to liftoff along the direction of leg spring, the value of which changes according to the real-time leg length.
Extension method. The 3D-HFC method mentioned above is developed assuming a 3D-SLIP equivalent model. However, ultimately, the control is intended for a quadruped robot. Then we are able to extend 3D-HFC approach as follows. First the control strategy for each joint could be chosen and coordinated in a way that achieves the equivalent virtual spring characteristics. Secondly, the multi-joint one-legged model can be expanded to the four-legged model. The control method of each leg is the same as the one-legged model. Then, a variety of gaits can be achieved by alternating motion of each leg. Figure 7 is the extension process of the model. As the Fig. 7b shown, the multi-joint one-legged model is composed of a mass, three connecting rods and a toe. The one-legged system contains four active control joints, namely the side-swing joint, hip joint, knee joint and ankle joint respectively. In the 3D-SLIP model, only side-swing joint and hip joint require active control. However, knee joint and ankle joint require additional active control in the one-legged model. Four joints should be cooperatively controlled to realize the equivalent virtual spring features. Figure 8 shows the extended schematic diagram of the multi-joint one-legged model. Joint angles and driving torque are defined as (θ 1 , θ 2 , θ 3 , θ 4 ) and (τ 1 , τ 2 , τ 3 , τ 4 ) respectively. Points O and P denote the CoM of the body and the foot endpoint. The idea of equivalent extension is to add a virtual spring between the point O and point P. The driving torques of the side-swing joint and hip joint are used to adjust the body posture during the stance phase and control the touchdown angle during the flight phase. The equivalent spring force and compensation force are achieved by torque control of the knee joint. Connecting rod BP is limited to the virtual spring direction through torque control of the ankle joint. The green dashed line between points O and P in the Fig. 8b represents the equivalent virtual leg spring. Its equivalent length, virtual stiffness and viscous damping are defined as (r S , k S , c S ).
Equivalent extension process is divided into flight phase and stance phase. In the flight phase, the side-swing angle ϕ TD and forward-swing angle θ TD of the virtual spring correspond to the angles in the side-swing plane and forward-swing plane of the 3D-SLIP model respectively. Theoretical values of these angles can be calculated by the 3D-HFC method in real time from the Eq. (5). The equivalent length of the virtual spring need to keep  www.nature.com/scientificreports/ the original length in the flight phase. The expressions of each joint angle can be obtained through the geometric relationship as follows: where r S0 represents the original length.
In the stance phase, the connecting rod BP points to the CoM along the direction of the virtual spring. Then the equivalent length r S can be described as follows: Similar to the 3D-SLIP model, equivalent force of the one-legged model also consists of the virtual spring force F S and the compensational force F E . F E can be calculated using the Eq. (11) with the 3D-HFC method. Considering the equivalent stiffness and damping, F S can be described as where v S represents the CoM velocity along the direction of the virtual spring. Meanwhile, the body attitude can be adjusted by the side-swing joint and hip joint. The torque expressions are the same as Eq. (6) following the body attitude control strategy of the 3D-HFC method. Then the expressions of each joint torque can be described as follows: where d S represents the moment arm between point A and the virtual spring OP, and θ 4,des represents the desired angle which can ensure the rod BP along the direction of the virtual spring. Both of them can be obtained through the geometric relationship as follows: Equations (12)-(17) provide the extension method from the equivalent model to the multi-joint one-legged model. Further, we can expand the control method to the four-legged model. In the case of the trot gait, four legs of the system can be divided into two groups, and each group includes two legs in a diagonal direction. The motion sequences of two legs in a group are approximately the same. Two groups of legs are usually in a different motion state, one is in the flight phase, and the other one is in the stance phase. Equations (12) and (15) respectively represent the control algorithm in the flight phase and stance phase.
If system body receives the impact disturbance, because of the same motion sequences, one group legs in the flight phase can be seen as an equivalent spring and adjust the touchdown point within the 3D-SLIP algorithm. And another group legs in the stance phase also can be seen as an equivalent spring and provide the compensation force and attitude control torques. Then the extended control methods can be implemented on the four-legged model.

Results
Simulation. The ADAMS-Simulink co-simulation technique is employed to verify the validity of the 3D-HFC strategy for 3D-SLIP model under external impact force. There are five different impact situations in Co-simulations: lateral impact when the SLIP system jumps vertically, lateral impact when the SLIP system runs forward, impact along with forwarding motion, compound impact in the running, and lateral impact for extension models. To meet the requirements of the real-time control system, the communication interval between ADAMS and Simulink in all simulations is 0.1 ms 44 . www.nature.com/scientificreports/ Modeling for simulation. The 3D-SLIP running model built in ADAMS is shown in Fig. 9. It includes the body, link 1, link 2, a spring with damper, and a toe. The coordinate systems and angles were defined identically to Fig. 3. The revolute joint 1 was placed between the body and link 1, and the revolute joint 2 was placed between link 1 and link 2. Between link 2 and leg toe, a translational joint was set to limit the movement direction of the leg spring. The compensational force was applied on the body CoM along the spring direction. The contact was defined between the leg toe and the ground. Moreover, a contact sensor was set to judge the touchdown event and liftoff event of the system. Except for the contact sensor, the sensing system also includes the body velocity sensor, spring length sensor, swing angle sensor, etc. The 3D-SLIP model parameters were referenced to the physiological index of physical dogs 47 , which are shown in Table 1. The Adams-Simulink co-simulation framework based on 3D-HFC is shown in Fig. 10.
Lateral impact when jumping. Vertically Simulation 1 (S1) was designed for the lateral impact of the SLIP system during vertical jumps. The initial apex height was set to 840 mm, while the initial forward velocity (along the x-axis) and lateral velocity (along the z-axis) were set to zero. An instantaneous lateral impact was applied on the body CoM along the z-axis at 0.92 s. The impulse value was about 20 kgm/s. The SLIP system did not fall down under the lateral impact. The initial state, including apex height, forward velocity and lateral velocity, was restored using the 3D-HFC method. The motion sequence diagram of this simulation is shown in Fig. 11. The red curve and the blue curve represent the motion trajectories of the body CoM and the toe in the lateral plane, respectively. Simulation results are shown in Fig. 12. The whole process of the SLIP system achieved the expected goal with the control of the 3D-HFC method according to the height curve (Fig. 12a). Although there is a little fluctuation after impact, the CoM apex height converges to the expected value very quickly. In contrast, the lateral velocity is more obviously affected after the impact. Even though the lateral velocity changed instantaneously, it can still converge rapidly to the desired level in two cycles. Figure 12b is the phase portrait, which shows that the return map can converge to closed curves, indicating the stable periodical motion achieved by the system. It means the system has been controlled at the target state, achieving the desired periodical motion. Impact in running. Simulation 2-4 (S2-S4) was designed for the lateral impact when SLIP system is running forward. In these simulations, the initial forward velocity of the model was set to 2.8 m/s to conduct a running  www.nature.com/scientificreports/ simulation on the flat terrain. The initial apex height and lateral velocity were set to 840 mm and zero, respectively. In S2, an instantaneous lateral impact was applied on the body CoM along the z-axis at 0.92 s, with an impulse value of about 20 kgm/s. In S3, by contrast, the direction of the impact force was changed from the z-axis to the x-axis. Considering the initial forward velocity, the impulse value was reduced to 10 kgm/s. In addition, these above two impacts were simultaneously applied on the body CoM at 0.92 s in S4. The impacts can be seen as different components of a spatial impact, respectively, on the x-axis and z-axis. The value and direction of the spatial impact can be defined by changing the values of the lateral impact and forward impact. The schematic diagram of different impacts is shown in Fig. 13. The locomotion of SLIP system is extended from 2-D plane to 3-D space as there is a lateral impact. The revolute joint and must be adjusted simultaneously to realize the dynamic stability control in 3-D space. The motion sequence diagram of S4 is shown in Fig. 13. Similarly, the red and blue curves represent the motion trajectories of the body CoM and the toe, respectively. On the X-Z projection plane, it can be clearly seen that the SLIP system deviates from the predetermined direction after impact, but soon recovers stability and continues to move forward along the predetermined direction.
The simulation results are shown in Fig. 14. Figure 14a-c show the curves of body CoM trajectory, forward and lateral velocities in S2-S4. It can be seen that the apex height is affected obviously by the impact along forward  www.nature.com/scientificreports/ motion, rather than the lateral impact. And the forward and lateral velocities are affected immediately by forward and lateral impacts respectively. Figure 14d-f show the apex height error and the forward velocity error in each step. When there is only a lateral impact, the errors of apex height and forward velocity are controlled within 5%. And when there is an impact alone the forward motion, apex height error and the forward velocity error can reach about 22% and 15% respectively under the unknown impact. However, the SLIP system will not fall down and can converge quickly to the expected state with an error less than 1% using 3D-HFC method. Figure 14g-i are the phase portraits and it can clearly find that there is a cycle obviously leaving the convergence range after impact. And after few steps, the return map is able to converge to closed curves, indicating the stable periodical motion achieved by the system. S4 is the most complex situation in these simulations, and the SLIP system can still be controlled at the target state and achieve dynamic stability motion. This means that the 3D-HFC control strategy has certain robustness and automatic adaptability, and is able to achieve the control target of the system's convergence to stability even under a spatial impact.
Different impact forces. In order to establish the performance limits of the control method, different impact forces were added on the body CoM by scanning. The value of the impulse started with 10 kgm/s, and increased by 10 kgm/s every time. Figure 15 shows the motion paths in Z-X plane under different impact forces. Due to different values of the lateral impact, the lateral velocities after impact are also different. When the impulse is set to 10 kgm/s, the maximum lateral displacement is 160 mm. The system can restore stability within one cycle after the impact and the final offset is 140 mm. When the impulse is set to 30 kgm/s, the system takes two cycles to restore stability, and the final offset increases to 250 mm. However, when the impulse is set to 50 kgm/s, the impact is too large and causes the system to fall after two cycles. This is beyond the capability of the 3D-HFC strategy. After several additional attempts, we found that the stability threshold of the impact was between 44 and 45 kgm/s. Figure 16a-c are the phase portraits under three different impact forces. With the increase of the impulse value, the convergence speed of the system reduces gradually until the system is in a state of divergence. We can more clearly see the effect of the control method from the phase portraits.   www.nature.com/scientificreports/ Lateral impact for extension models. In order to verify the effectiveness of the extension methods, the co-simulations were also performed on the one-legged model and four-legged model. Initial conditions and control targets were the same as the S2. Similarly, a uniform instantaneous lateral impact was added on the body CoM of the 3D-SLIP model, one-legged model and four-legged model respectively. The value of the impulse was about 20 kgm/s. Figure 17a shows the motion sequence diagrams for the one-legged model. In the XY projection plane, the one-legged model can achieve the periodic forward motion similar to the 3D-SLIP model. And in the XZ projection plane, we can also see that the one-legged system deviates from the intended direction after impact, but soon restores stability and continues to move forward along the intended direction. Figure 17b shows the motion sequence diagrams for the four-legged model. The four-legged system can also move forward steadily and periodically within the extension method. Figure 18 shows the simulation results of the comparison for different models. Obviously, the characteristics of the CoM trajectory and lateral velocity between three models are quite similar under the lateral impact. This means that the extension method has certain effectiveness for the one-legged model and the four-legged model.

Prototype experiment.
A further goal is to demonstrate the 3D-HFC method on a real legged robot prototype. Although there are certain differences between the equivalent model and the real robot system, the lateral impact experiment is designed for a real robot system to demonstrate and further improve the effectiveness of the balance control method through experimental results.
In order to enhance the performance of quadruped robot, a prototype named MBBOT was developed by a team including researchers from Harbin Institute of Technology (HIT), Huazhong University of Science and Technology (HUST), Shenyang Institute of Automation (SIA), Chinese Academy of Sciences, University of Electronic Science and Technology of China (UESTC) and the Academy of Armored Forces Engineering (AAFE).    59 and so on. After studying these robot systems, we propose our design for the quadruped platform. As shown in Fig. 19a, the MBBOT consists of a torso module and four leg modules. The torso module is used to fix the power system and central control unit, and to connect the four leg modules as a platform. Each leg module contains four active rotary joints and a passive spring, and has four active DOFs and a passive DOF. The four active rotary joints correspond to the Side-swing Joint, Hip Joint, Knee Joint and Ankle Joint respectively, and they are driven by identical hydraulic cylinders. All 16 joint-actuators are powered by the external hydraulic pump which is also fixed on the torso module. There are 41 sensors on the robot prototype, such as displacement sensors and load cells on the hydraulic cylinders and springs, 3-component force sensors in feet and an Inertial Measurement Units (IMU) on torso. The controller of the MBBOT is the PC-3363 bus based on PC/104. It is a small size, low power consumption, superior performance single board computer (SBC) with built-in floating point operation unit (FPU). PC-3363 board adopts the stacking structure, so it is very convenient to expand function modules. The details of the controller are shown in Table 2. The total weight of the MBBOT is about 150 kg, and the mass ratio between torso and single leg is about 10:1.
The lateral impact experiment platform is shown in Fig. 19b. A sandbag weighing about 40 kg is hung on the bracing frame. The lateral impact force is generated by the sandbag as a swinging pendulum. The sandbag impacts the MBBOT at the lowest swing point and the impact value can be estimated by the relative vertical distance between the initial swing point and the lowest swing point. For example, the relative vertical distance is about 0.11 m in the lateral impact experiment, and the instantaneous impulse value can reach 60 kgm/s. Protection frame is used to prevent the robot falling down if the control algorithm fails during the lateral impact experiment. The action sequence diagram of the impact experiment is shown in Fig. 20.  www.nature.com/scientificreports/ The experimental results are shown in Fig. 21. The initial state of the MBBOT is stepping there waiting for the impact from the swinging sandbag. A moment later, the sandbag swings to the lowest point and impacts on the MBBOT. The IMU is detecting the body attitude angles and the change of the acceleration, and then the planning joint angles can be calculated by the extension method based on the 3D-HFC algorithm. Figure 21a is the joint angles of the back-left leg (BL-leg) during the impact progress. Solid lines and dotted lines represent the measured values (MV) and the desired value (DV). Figure 21b shows the vertical force at all legs foot-end for each duration. During the stance phase, the vertical force is high. And in contrast, vertical force is low during the flight phase. In the graph, the gait-pattern shows obvious symmetric characteristics. The attitude angles are shown in the Fig. 21c. We can see that the attitude angle γ is most affected by the lateral impact and the maximum deflection angle can reach 11 • . And Fig. 21d shows the range of the attitude angle for the experiment and simulation. By the analysis of the results, attitude angles in simulation oscillated around the desired angle and had a smaller fluctuation range. The maximum angular deviation was about 2.8 • . In contrast, attitude angles in experiment more deviated from the desired angle, and had a greater fluctuation range. The maximum angular deviation was about 11.1 • . The differences of control effect may be caused by the leg quality of the prototype robot. Another possible reason was that the sandbag didn't precisely impact on the CoM of the prototype. However, the MBBOT could restore stability after the impact with the extension method. This experiment verifies the feasibility and effectiveness of the extension based on 3D-HFC algorithm under the instantaneous lateral impact.

Discussion
A dynamic stability control method is proposed for SLIP running systems in 3D space in this study. In case of an unknown instantaneous force impact, the SLIP system is able to quickly restore stability without falling from impact disturbance. As a result of hybrid feedback control, the 3D-HFC method has a strong ability of fault tolerance control and environmental adaptation performance. The energy loss in each step is taken into consideration in the controller to improve the stability. The co-simulations verify the good control effect of system under lateral impact when jumping vertically, lateral impact when running forward, impact along forward motion and www.nature.com/scientificreports/ compound impact in running. This method has also been extended to apply on the real legged robot platform, and a lateral impact experiment has been designed and implemented.
Since the ideal equivalent SLIP model ignores the mass of the leg, the spatial rotational kinetic energy is not taken into account in the theoretical dynamic model. However, the leg mass is not zero in the real system, and the body attitude angle will change with the motion of the leg. This is the reason why the 3D-HFC method increases the body attitude control module. Fortunately, the mass ratio between body and single leg is about 10:1, the influence of the mass and inertia of the leg can be reduced or eliminated by the body attitude adjustment control. And the extension method also could be used in the real robot prototype. But in the further work, the improved dynamic model will be intended to establish considering the mass and inertia of the leg.